function mpr = getMPR_macro_spec(output,theta2,AVarTheta,nx,switch_h0,switch_hx)


% Version 1: MPR = (lambda0 + lambdax)
% Version 2: MPR = (lambda0 + lambdax)*inv(sigma)
version = 2;

[h0_1,h0_2,hx_1,hx_2,sigma] = unfoldTheta2_threshold_macro(theta2,nx);
if version == 1
    lambda0_1 = sigma\(h0_1-output.model.phi*output.model.mu);
    lambda0_2 = sigma\(h0_2-output.model.phi*output.model.mu);
    lambdax_1 = sigma\(hx_1-(eye(nx)-output.model.phi));
    lambdax_2 = sigma\(hx_2-(eye(nx)-output.model.phi));
elseif version == 2
    lambda0_1 = h0_1-output.model.phi*output.model.mu;
    lambda0_2 = h0_2-output.model.phi*output.model.mu;
    lambdax_1 = (hx_1-(eye(nx)-output.model.phi));
    lambdax_2 = (hx_2-(eye(nx)-output.model.phi));
end

if version == 1
    % Computing the standard errors by simulation (ignoring uncertainty about
    % theta1, which is ok for T going to infinitey because then N*T goes much
    % faster to infinity given N also goes to infinity)
    numSim = 50000;
    rng(1,'twister')
    lambda0_1Draw = zeros(nx,numSim);
    lambda0_2Draw = zeros(nx,numSim);
    lambdax_1Draw = zeros(nx,nx,numSim);
    lambdax_2Draw = zeros(nx,nx,numSim);
    namesTheta2Reduced = AVarTheta.names;
    numTheta2Reduced   = length(namesTheta2Reduced);
    for i=1:numTheta2Reduced
        theta2Reduced.(namesTheta2Reduced{i}) = theta2.(namesTheta2Reduced{i});
    end
    [Stheta2Reduced,check] = chol(AVarTheta.VarRobust,'lower');
    if check ~= 0
        mpr.lambda0_1       = lambda0_1;
        mpr.lambda0_2       = lambda0_2;
        mpr.lambdax_1       = lambdax_1;
        mpr.lambdax_2       = lambdax_2;
        mpr.SE = NaN;
        return
    end
    for s=1:numSim
        theta2ReducedDrawValues = struct2array(theta2Reduced)'+Stheta2Reduced*randn(numTheta2Reduced,1);
        theta2ReducedDraw       = values2struct(theta2ReducedDrawValues,namesTheta2Reduced);
        % Dublicating parameters which are common across the two regimes
        if any(switch_h0 == 0)
            for i=1:nx
                if switch_h0(i,1) == 0
                    theta2ReducedDraw.(['h0',num2str(i),'_2']) = theta2ReducedDraw.(['h0',num2str(i),'_1']);
                end
            end
        end
        if any(any(switch_hx == 0))
            for i=1:nx
                for j=1:nx
                    if switch_hx(i,j) == 0
                        theta2ReducedDraw.(['hx',num2str(i),num2str(j),'_2']) = theta2ReducedDraw.(['hx',num2str(i),num2str(j),'_1']);
                    end
                end
            end
        end
        [h0_1,h0_2,hx_1,hx_2,sigma] = unfoldTheta2_threshold(theta2ReducedDraw,nx);
        % Intercepts for MPR
        %lambda0_1Draw(:,s) = sigma\h0_1;
        %lambda0_2Draw(:,s) = sigma\h0_2;
        lambda0_1Draw(:,s) = sigma\(h0_1-output.model.phi*output.model.mu);
        lambda0_2Draw(:,s) = sigma\(h0_2-output.model.phi*output.model.mu);
        % Loadings for MPR
        %lambdax_1Draw(:,:,s) = sigma\(hx_1-(eye(nx)-outputStep1.model.phi));
        %lambdax_2Draw(:,:,s) = sigma\(hx_2-(eye(nx)-outputStep1.model.phi));
        lambdax_1Draw(:,:,s) = sigma\(hx_1-(eye(nx)-output.model.phi));
        lambdax_2Draw(:,:,s) = sigma\(hx_2-(eye(nx)-output.model.phi));
    end
    lambda0_1_SE = std(lambda0_1Draw,0,2);
    lambda0_2_SE = std(lambda0_2Draw,0,2);
    lambdax_1_SE = std(lambdax_1Draw,0,3);
    lambdax_2_SE = std(lambdax_2Draw,0,3);
elseif version == 2
    theta2se = struct();
    for i=1:length(AVarTheta.names)
        name = AVarTheta.names{i};
        theta2se.(name) = sqrt(AVarTheta.VarRobust(i,i));
    end
    % Dublicating parameters which are common across the two regimes
    if any(switch_h0 == 0)
        for i=1:nx
            if switch_h0(i,1) == 0
                theta2se.(['h0',num2str(i),'_2']) = theta2se.(['h0',num2str(i),'_1']);
            end
        end
    end
    if any(any(switch_hx == 0))
        for i=1:nx
            for j=1:nx
                if switch_hx(i,j) == 0
                    theta2se.(['hx',num2str(i),num2str(j),'_2']) = theta2se.(['hx',num2str(i),num2str(j),'_1']);
                end
            end
        end
    end
    [h0_1se,h0_2se,hx_1se,hx_2se] = unfoldTheta2_threshold_macro(theta2se,nx);
    lambda0_1_SE = h0_1se;
    lambda0_2_SE = h0_2se;
    lambdax_1_SE = hx_1se;
    lambdax_2_SE = hx_2se;   
end

% Computing t-stats
lambda0_1_tstat = lambda0_1./lambda0_1_SE;
lambda0_2_tstat = lambda0_2./lambda0_2_SE;
lambdax_1_tstat = lambdax_1./lambdax_1_SE;
lambdax_2_tstat = lambdax_2./lambdax_2_SE;

%% saving results
mpr.lambda0_1       = lambda0_1;
mpr.lambda0_2       = lambda0_2;
mpr.lambdax_1       = lambdax_1;
mpr.lambdax_2       = lambdax_2;
mpr.lambda0_1_SE    = lambda0_1_SE;
mpr.lambda0_2_SE    = lambda0_2_SE;
mpr.lambdax_1_SE    = lambdax_1_SE;
mpr.lambdax_2_SE    = lambdax_2_SE;
mpr.lambda0_1_tstat = lambda0_1_tstat;
mpr.lambda0_2_tstat = lambda0_2_tstat;
mpr.lambdax_1_tstat = lambdax_1_tstat;
mpr.lambdax_2_tstat = lambdax_2_tstat;
end